% Simulate surprise closure of Newberry hunt. Assume those awarded a permit
% to a Newberry hunt do not go, but no one else is affected.

clear

rng(1196)

% Load data, assign parameters
load('sigma0save.mat')
pit = [thetastar(4) thetastar(5) 1-thetastar(4)-thetastar(5)];              % Type probabilities
s = thetastar(end);                                                         % Cognitive cost parameter
mu = thetastar(1);                                                          % Travel cost parameter
eta = [thetastar(2) thetastar(3) 0];                                        % Type fixed effects
chi = [0 thetastar(6:end-2) 0];                                             % Hunt baseline utilities
phi0(end,:) = zeros(1,K);
tc = [zeros(size(tc,1),1) tc zeros(size(tc,1),1)];

% Specify type, travel cost distributions from model estimates and data
typedist = makedist('Multinomial',pit);                                         
tcdist = makedist('Multinomial',PTC);

% Set sites to be closed (Newberry hunts)
closed = 17:19; 

% Initialize result vectors, set number of iterations
iter = 1e6;
tctracker = zeros(iter,1);
typetracker = zeros(iter,1);
pptracker = zeros(iter,1);
jstartracker = zeros(iter,1);
phiv0 = zeros(iter,1);
phiv1 = zeros(iter,1);

% Run simulation
parfor idx1 = 1:iter
    
    % Draw travel cost profile, type
    tcdraw = random(tcdist);
    tctracker(idx1) = tcdraw;
    typedraw = random(typedist);
    typetracker(idx1) = typedraw;
    
    % Define pref pt distribution conditional on type, tc profile, then
    % draw pref pt stock
    ppdist = makedist('Multinomial',zeta0(:,:,tcdraw,typedraw)');
    ppdraw = random(ppdist);
    pptracker(idx1) = ppdraw;
    
    % Define type-1 EV and exponential error terms (from hunting (e) and cognitive cost (w))
    e = -log(-log(rand(H-1,1)));
    w = log(1 - rand())/-s;

    error = phi0(:,ppdraw).*[0; e; 0] + [zeros(H,1); w];
    [~,jstar] = max(Vb0(:,ppdraw,tcdraw,typedraw) + error);
    jstartracker(idx1) = jstar;
    dropped = sum(jstar==closed) > 0;
    
    % Calculate flow utility before, after surprise closure
    phiv0(idx1) = phi0(jstar,ppdraw)*(chi(jstar) + eta(typedraw) ...
        + mu*tc(tcdraw,jstar)) + error(jstar);
    phiv1(idx1) = (1 - dropped)*phiv0(idx1); 
        
end

WTP_Newberry = mean(phiv1(jstartracker==17|jstartracker==18|jstartracker==19) ...
    - phiv0(jstartracker==17|jstartracker==18|jstartracker==19))/abs(mu/200)
% mean(phiv0)/abs(mu/200)

WTP_Not_Newberry = mean(phiv1(jstartracker~=17&jstartracker~=18&jstartracker~=19) ...
    - phiv0(jstartracker~=17&jstartracker~=18&jstartracker~=19))/abs(mu/200)

Newberry_prop = sum(jstartracker==17|jstartracker==18|jstartracker==19)...
    /size(phiv0,1)

save('Results.mat')
